The Asymptotic Falloff of Local Waveform Measurements in Numerical Relativity 
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We examine current numerical relativity computations of gravitational waves, which typically determine the 
asymptotic waves at infinity by extrapolation from finite (small) radii. Using simulations of a black hole binary 
with accurate wave extraction at r = 1000M, we show that extrapolations from the near-zone are self-consistent 
in approximating measurements at this radius, although with a somewhat reduced accuracy. We verify that ^4 is 
the dominant asymptotic contribution to the gravitational energy (as required by the peeling theorem) but point 
out that gauge effects may complicate the interpretation of the other Weyl components. 
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/. Introduction. - Numerical relativity has made great strides 
in recent years in the solution of the binary black hole (BH) 
problem. Since the original breakthroughs by Pretorius [ 1 1 
and the moving puncture approach J2] 0, the calculation of 
long, accurate gravitational waveforms (GWs) has become an 
almost routine procedure. It is particularly satisfying that a 
variety of methods (numerical methods, formulations of the 
Einstein equations, wave extraction techniques) are in use and 
have been shown to produce consistent results (eg. [4|). 

Certain systematic errors, however, are difficult to estimate. 
In particular, current methods measure GWs at finite radii and 
extrapolate the results to infinity. This extrapolation has been 
identified as one of the largest remaining sources of systematic 
error within current extraction techniques, particularly dur- 
ing the merger and ring-down (5j |6). Potential ambiguities 
arise particularly at small radii where gauge as well as nonlin- 
ear near-zone effects may dominate the expected polynomial 
falloff of the amplitude. 

In this paper, we verify the extrapolation procedure for 
GWs by performing accurate wave extractions at large radii, 
out to r = 1000M from the source (where M is the mass 
of the spacetime). The waveforms have been calculated using 
a new hybrid multi-patch/mesh-refinement algorithm, which 
allows for an efficient discretisation of the wave zone so that 
high accuracy can be obtained to large radii. We find that 
the measured waves between r — 100M and r — lOOOAf 
are convergent and of good enough quality to extrapolate the 
phase and amplitude accurately by low-order polynomial ex- 
pansions. The measurements at r = 1000M can be estimated 
to within 0.04% in amplitude and O.OOlrad in phase, if the 
measurements out to r = 600Af are used in the extrapolation. 
This is true over the course of the evolution, including 8 orbits 
of inspiral, the merger, and ring-down. If only measurements 
within r = 220M are used, as is common, then the errors 
increase by an order of magnitude. 

Finally, we note that the gravitational radiation is normally 
associated with the leading order term in the falloff of the 
spacetime curvature. By the peeling theorem, we expect this 
to be a polynomial in 1/r whose leading coefficient is the 
Weyl component ip 4 . By measuring all of the Weyl curvature 



components, we have been able to establish their respective 
falloff rates, and verify that ip4 is indeed the leading order co- 
efficient with the expected 1/r falloff rate. The exponents for 
the ipi and ipo components are less clear, however, and suggest 
that local gauge effects influence their computation. 
//. Computational setup. - A key feature of the calculations 
performed here is the accuracy which we are able to achieve at 
large radii from the source through the use of a newly imple- 
mented numerical scheme. The code makes use of finite dif- 
ferences and standard mesh-refinement techniques, but incor- 
porates the use of multiple grid patches to cover the spacetime 
with flexible adapted local coordinates. For the binary BH in- 
spiral considered here, we consider two regions, depicted in 
Fig[T] In the near-zone region where the BHs orbit, we discre- 
tise the spacetime using standard Cartesian grids, applying 2: 1 
Berger-Oliger mesh refinement in order to increase the resolu- 
tion around each body |7|. In the wave-zone, however, the dy- 
namical fields are essentially radially propagating waves. We 
cover this zone with six overlapping patches, each of which 
incorporate a local radial coordinate r and transverse angular 
coordinates (p, a). The use of six patches avoids the problem 
of a coordinate singularity on the axis of a single spherical- 
polar coordinate system, as well as providing a more uniform 
angular resolution over the sphere. The particular coordinates 
which we have implemented are the "inflated cube" coordi- 
nates, given explicitly in [8 1. 

Derivatives on each grid are locally computed using stan- 
dard finite differences at 8th-order. Data is passed between 
patches by interpolation, typically via centred 5th-order La- 
grange polynomials. Each patch is surrounded by a boundary 
zone which is populated with data mapped from the neigh- 
bouring patch, so that derivatives can be calculated up to the 
patch edge without the need for off-centred stencils. A 4th- 
order Runge-Kutta integrator is used to evolve the solution. 

We write the Einstein equations in the commonly used 
BSSNOK form (H US EE E), adopting the particular vari- 
ation proposed by 1131 . whereby the usual variable <fi — 
log 7/I2 is replaced by w = 7~ 2 (with 7 the 3-metric deter- 
minant). Gauges are the commonly used 1 + log and T-driver 
conditions with advection terms (2] [3] ED . 
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FIG. 1: Schematic depiction of the grid structure in the z — plane. 
Four radial grids surround the equator. The right shows an expanded 
inset of the Cartesian grid (shaded) covering the near-zone around 
the individual BHs. 



The GWs are measured by evaluating the Weyl curvature 
tensor components, Cap^s, in a null frame (1, n, m, m), ori- 
ented so that the outgoing vector, 1, points along the coordi- 
nate f direction while the other vectors determine an orthonor- 
mal null frame in the local metric. The independent curvature 
components ipo . . . -04 are determined by the standard projec- 
tions lfl5l . The Weyl components are evaluated on spheres of 
fixed coordinate radius and projected onto a basis of spherical 
harmonics, -nXi m - For binary black hole situation considered 
here, the dominant mode is I = 2, m = 2, which is used for 
the results presented here. 

To establish the numerical accuracy, we have performed 
evolutions at three different resolutions, and find that the re- 
sults converge at 4th-order during the merger, and close to 
8th-order during the inspiral, in both amplitude and phase. (A 
more complete description of the code implementation and 
tests will be given elsewhere |16j].) Here we present results 
based on the highest resolution evolved, for which the spa- 
tial resolution for all of the GW measurements is uniformly 
h r = 0.64A7 in the radial direction, and h± ~ 3° in the an- 
gular directions. GW measurements are taken every 0.144M. 
III. Extrapolation of waveforms. - We have evolved an equal- 
mass, non-spinning binary from separation d/M = 11.0 
through approximately 8 orbits (a physical time of around 
1360Af), merger and ring-down. The masses of the punctures 
are set to m — (0.4872) and are initially placed on the a>axis 
with momenta p = (±0.0903, TO. 000728, 0), giving the ini- 
tial slice an ADM mass M = 0.99052. These initial data pa- 
rameters were determined using a post-Newtonian evolution 
from large initial separation, following the procedure outlined 
in ifTTI . with the conservative part of the Hamiltonian accurate 
to 3PN and radiation-reaction to 3.5PN, and determine orbits 
with eccentricity less than 0.2%. 

We measure the Weyl components ipo . . .ipi every 20M 
from r = 100 A/ to 300 M, then at 400M, 500 Af, 600M and 
lOOOAf. The radial grid structure in the wave zone allows us 
to extend the outer boundary of the grid at relatively little cost 
compared to Cartesian codes. For the runs presented here, it is 
placed at r = 3600A/ with a resolution of dr = 2.56Af at the 



outer boundary so that the I = 2, m = 2 mode is reasonably 
well resolved throughout the grid. This allows for 2600A/ of 
evolution time before a physical or constraint violating mode 
traveling at the speed of light can reach the outermost detec- 
tor at r = 1000M. That is, the outer boundary is effectively 
causally disconnected from the wave measurements presented 
in this paper. 

The Weyl components ipj = Aje 1 ^' are assumed to fall-off 
as a function of radius according to 

Mr,n = E Mr,n = E 

i=0 i=0 

(1) 

The r coordinate is that of the simulation coordinates, which 
we find to differ by at most 0.1% from the areal radius. The 
GWs are expressed in terms of the retarded time t* = t — r* 
where t is the coordinate time and 7-* = r + 2M ln[r/(2Af) — 
1] is the tortoise coordinate [5]. We do not offset the retarded 
time to align the peaks of the waveforms. 

It is generally difficult to estimate the error incurred when 
extrapolating. Given the data at r = 1000M, we can at- 
tempt to gauge an optimal choice of extrapolation parame- 
ters by attempting to estimate this data from the measure- 
ments at smaller radii. As test cases, we construct extrapola- 
tions using four different sets of radii, e x S {100A/, 200A/}, 
e 2 e {160M,280M}, e 3 G {200Af,300Af} and e 4 € 
{260A/, 600A/}. Each of these incorporates 6 data points, 
which over-determines low-order polynomials. We evaluate 
the extrapolation coefficients by a least-squares fit to these 
points, which can be important in removing spurious oscil- 
lations that may arise fitting high-order polynomials to noisy 
data. 

Fitting the amplitude using various polynomial orders, tia, 
suggests that in all cases = 3 is optimal in predicting the 
amplitude of the measured wave at r — 1000 M, with an er- 
ror of approximately 0.02% for e 4 and 0.2% for e\. For the 
phase, we find that = 3 minimises the error, at 6 x 10 _4 rad 
and 5 x 10~ 3 rad for e 4 and ei, respectively. In both amplitude 
and phase, we note that the error is reduced significantly if the 
outermost data e 4 is used. In Fig.|2]we display the error in es- 
timating the r = lOOOAf data using each of the extrapolations 
at the optimal order. The maximum errors in both amplitude 
and phase tend to occur during the late inspiral (t* = — 200 Af 
to t* — 0M) and ring-down, although there is no sign of a 
rapid growth of error during this phase. 

The corresponding extrapolations to r — > oo shows very 
similar behaviour. In Fig. [3] we have compared each of the 
extrapolations with an extrapolation obtained by including the 
r = 1000M data (e 5 € {280Af, 1000A/}), evaluated at r -> 
oo. The outermost extrapolations differ by at most A A = 
0.03% and 5<j> — 0.003rad over the course of the evolution. 
IV. Peeling properties. - The interpretation of ip4 as the ra- 
diated gravitational energy is a consequence of the "peeling" 
property, which states that for asymptotically flat spacetimes 
at large radii, the Weyl curvature tensor can be represented 
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FIG. 2: Error in the extrapolated amplitude (top panel) and phase 
(bottom panel) of the I — 2, m = 2 component of 1p4d.tr = 1000A/ 
as computed by extrapolations ei . . . e4 (defined in the text). 
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FIG. 3: Differences between the r < 600M extrapolations with an 
extrapolation including r — 1000M data (es) evaluated in the limit 



schematically as 
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— + -o + -3+ -z + -ir + 0(l/r°). (2) 



That is, each component of the Weyl tensor falls off at a 
known fixed rate, and at large radii, ipi is the dominant com- 
ponent. At future null infinity, J + , it can be related to M, 
the change in energy of the spacetime. We note, however, that 
the peeling theorem involves a number of restrictions on the 
asymptotic form of the spacetime, and the coordinates which 
are used there. A rigorous connection between finite radius 
measurements and the asymptotic properties of the spacetime 
at J + is difficult to make. 

Given the importance of the falloff of the curvature in the 
identification of ^4 with the GW, it is useful to examine the 



behaviour of the other Weyl components measured by the sim- 
ulation. In Fig. [4j we have plotted their falloff as a function 
of coordinate radius. The time series data for each component 
is mapped to a scalar by integrating the amplitude over the in- 
terval t e [-800M, 50M]. (Alternatively, one could obtain a 
scalar by taking the measurements at a point such as the wave- 
form peak. A similar plot results, but the averaging effect of 
the integral reduces local noise slightly.) 

For the cases of ^4 and ^3, we find that a straight line can 
be fitted to each of the components, indicating a consistent 
exponent, with measured values of —0.99 and —1.99 respec- 
tively, and a rather good agreement with Eq. Q. Due to its 
small amplitude, the tp2 measurement is dominated by numer- 
ical noise beyond a certain radius (clear from examination of 
the time-series data), and as a result, the curve veers from a 
straight line. However, if we fit a straight line to the five data 
points from r < 200A/, we find an exponent, —2.99, again 
agreeing well with the expectation. 

The ipi and ipo components present an interesting situation. 
Particularly in the case of ipo, the amplitude is large enough 
that a clear signal is present (of almost of the same ampli- 
tude as ^3). The falloff, however, is of order —2.00, rather 
than the —5 which the peeling theorem requires. Further, we 
note that the mode propagates outwards with a peak coinci- 
dent with that of 7/14, in contrast to the interpretation of ipo as 
an "ingoing" component of radiation. 

A possible explanation is that metric perturbations cause 
oscillations in the frame in which the components are mea- 
sured. As described above, we define the null frame only with 
reference to the local space and time coordinates. Attempts to 
modify the falloff of tpo via frame rotations (spin-boosts and 
null rotations) did not preserve the falloff of the other com- 
ponents. However, other gauge effects are likely present. We 
note that measurements are taken on spheres defined by the 
grid coordinates. The areal radius of these spheres exhibits 
small (on the order of 0.1%) oscillations in the I = 2, m — 2 
mode. The finite-radius ^4 measurement is known to be sus- 
ceptible to pure gauge effects such as the presence of a non- 
zero shift vector, which can produce spurious GW signals in 
static spacetimes ifTHIl . Though these effects are small, so are 
the values of tpi and tpo and thus correspondingly sensitive 
compared to that of the dominant component. 
V. Discussion. - We have demonstrated a number of features 
related to the measurement of GWs at finite radius. Our re- 
sults suggest that polynomial extrapolation of the ip4 com- 
ponent from small radii can provide an accurate model for 
estimating the measurements at larger radii. Data measured 
within r = 200M of the source have an error in amplitude 
and phase of AA ~ 0.2% and A(j> ~ 5 x 10~ 3 rad throughout 
the evolution (including merger and ring-down) compared to 
the measurement at r = 1000M . This provides an important 
check on numerical relativity measurements, which typically 
extrapolate from r < 200Af. Larger radius measurements 
do, however, improve the extrapolation, and errors can be re- 
duced by a further order of magnitude if data to r = 6Q0M 
is included. We also note that while ip4 is dominated by the 
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reducing boundary errors. 
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FIG. 4: The radial falloff of the Weyl components. Lines are linear 
least-squares fits to all of the points of lpo> ip3, an d t[>4, and the r < 
200M points for ipi and t/>2- Measured slopes are listed in the legend. 

1 jr term beyond r = 300Af , at smaller radii the higher order 
terms have a much larger contribution. 

While the falloff of ip 4 is the leading order contribution to 
the curvature as expected, the results of Fig.|4]suggest that the 
picture may be more complicated for the other Weyl compo- 
nents, and that care should be taken in interpreting local vari- 
ables according to their expected asymptotic properties. With 
respect to the peeling property, we note that the asymptotic 
evolution of generic initial data sets are likely to include poly- 
homogenous terms, and even the situation for Bowen-York 
data with linear momentum is not clear |[T9l l20ll . While it 
seems most likely that the effect we have observed in ipo and 
ipi is related to gauge, this provides a strong caution regard- 
ing physical predictions based on these quantities. Alternative 
gauge conditions may alleviate (or exaggerate) the issues we 
have noted. 

The mapping of finite radius results to asymptotic values at 
J + needs to be considered with some care. While our results 
suggest that the procedure of extrapolation is self-consistent 
and can be used to estimate the results that would be obtained 
by direct measurement at large radii, they do not establish 
the identification of the extrapolated quantities with quanti- 
ties that would be measured at J + . In a related paper Ell , 
we have demonstrated that the extrapolation procedure does 
in fact reproduce results obtained at J + to high accuracy, 
though a small systematic error does remain. A number of 
local corrections have been proposed to improve the rigor of 
■04 measurements at finite radius (cf. [22, 23] [24]). 

Finally, we note that our accurate measurements at r = 
lOOOAf are a result of a new computational infrastructure 
making use of adapted coordinate grids in conjunction with 
a finite difference, moving-puncture scheme. This is the first 
demonstration that such methods can produce stable evolu- 
tions for dynamical spacetimes. The efficiency gains allow 
the wave zone to be covered with sufficient resolution to very 
large radii (3600M in this case), which has been a crucial in 
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